function [nodes, weights] = legendre(lo, hi, num_nodes)

  
%  usage:  [nodes, weights] = legendre_quad(lo,hi,num_nodes)
%
%  Calculates the Legendre polynomial nodes and weights 
%  for given range and num_nodes.
%  From NR p. 152
  
  tol = (eps);
  m=(num_nodes+1)/2;
  xm=0.5*(hi+lo);
  xl=0.5*(hi-lo);
  
  
  ros     = max(length(hi),length(lo));
  nodes   = zeros( ros, num_nodes );
  weights = nodes;
  for i=1:m
  
    z=cos(pi*(i-0.25)/(num_nodes+0.5)); 
     % Initial guess at z
    z1 = z+.1;
    while abs(z-z1) > tol                 
      % then refine Newton's method...
      p1=1.0;
      p2=0.0;
      
      for j=1:num_nodes
     	p3=p2;
	    p2=p1;
	    p1=((2*j-1)*z*p2-(j-1)*p3)/j;
      end
      
      pp=num_nodes*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
      
    end
      
    nodes(:, i) = xm-xl*z;
    nodes(:, num_nodes+1-i) = xm+xl*z;
    weights(:, i) = 2*xl/((1-z*z)*pp*pp);
    weights(:, num_nodes+1-i) = weights(:, i);
 
  end
  
 %if ros==1
      nodes   = nodes';
      weights = weights';
 %end
